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Abstract 

The  calculation  of  neutron  flux  distribution  and  growth 
rate  for  small,  spherically  symmetric  systems  usually  requires 
extensive  computing  time  on  the  largest  machines.  To  minimize 
computing  time,  a  compromise  between  the  simplicity  of 
diffusion  theory  and  the  accuracy  of  transport  theory  is 
needed.  The  Serber-Wilson  method,  Feynman's  method,  and  early 
flux  synthesis  methods  are  used  as  the  foundation  for  integral 
integral  equation  synthesis  (IES)  which  is  an  approximate, 
numerical  technique  for  obtaining  the  spatial  and  energy 
neutron  flux  distributions  in  multiplying  systems.  In  IES, 
the  integral  form  of  the  neutron  transport  equation  is 
specialized  to  spatial  dependence  only,  and  then  solved 
numerically  for  the  two  lowest  order  eignefunctions .  Similar 
specialization  to  energy  dependence  only  yields  a  second  set 
of  trial  eigenfunctions.  Using  standard  perturbation  methods, 
the  two  sets  of  trial  eigenfunctions  are  synthesized  into  a 
single,  two-dimensional  solution.  The  IES  technique  was 
used  to  calculate  the  flux  and  multiplication  factor,  k  , 
of  the  critical  plutonium  sphere,  Jezebel.  Results  for  k 
agreed  to  within  0.01%  of  published  values,  whereas  the 
spatial  flux,  when  normalized  at  the  center,  agreed  to  within 
8%  at  the  outer  assembly  boundary.  The  Jezebel  calculation 
using  IES  required  about  90  seconds  CPU  time  on  an  IBM  360/75. 
Highly  sophisticated  codes  require  approximately  ten  minutes 
of  CDC  7600  CPU  time  to  compute  the  Jezebel  flux  and  growth 
rate. 
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INTEGRAL  EQUATION  SPACE-ENERGY  FLUX 
SYNTHESIS  FOR  SPHERICAL  SYSTEMS 

I.  Introduction 

The  calculations  of  neutron  flux  distributions  and 
neutron  growth  rate  (a)  for  multiplying  systems  currently 
require  extensive  computing  time  on  the  largest  machines. 
Because  of  the  immense  amount  of  detail  in  the  dependence 
of  cross  sections  on  neutron  energy  for  fissile  nuclei, 
there  is  no  possibility  of  obtaining  exact  solutions  to  the 
energy-dependent  neutron  transport  equation  for  general 
problems  (Ref  2:48).  The  usual  approach  is  to  obtain  numer¬ 
ical  solutions  to  the  Boltzmann  transport  equation.  Even 
for  one -dimensional  systems,  this  integro -differential 
equation  expresses  the  neutron  flux  as  a  function  of  four 
independent  variables:  position,  energy,  direction,  and 
time.  Accurate  calculations  of  neutron  transport  are  often 
made  by  finite  differencing  the  first  three  of  these  variables. 

A  typical  neutron  transport  problem  might  involve  tens 
of  spatial  cells  and  perhaps  ten  energy  groups  and  ten 
directions.  Differencing  this  problem  then  results  in  a 
mesh  composed  of  approximately  a  thousand  "cells;"  each 
with  a  specified  position,  energy,  and  neutron  direction. 

The  calculation,  however,  is  entirely  nonlocal  in  the  sense 
that  neutrons  in  any  cell  can  influence  neutrons  in  any 
other  cell.  Thus,  to  compute  a  neutron  flux  value  in  one 
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cell  requires  the  calculation  of  interactions  between  that 
cell  and  the  other  999  cells.  For  a  typical  problem  of  a 
thousand  cells  then,  about  a  million  interactions  must  be 
calculated  in  order  to  specify  the  neutron  flux  for  a  given 
time.  In  a  time  dependent  problem,  these  calculations 
must  be  repeated  each  time  the  clock  is  advanced  one  step. 
Clearly,  the  usual  approach  to  numerically  solving  the 
Boltzmann  transport  equation  requires  extensive  computing 
time  in  addition  to  extensive  computer  storage  capability. 

If  all  spatial  dimensions  of  the  system  under  consider¬ 
ation  are  large  compared  to  the  applicable  neutron  mean  free 
paths,  it  is  possible  to  use  diffusion  theory  to  calculate 
the  neutron  flux  and  growth  rate  for  the  system.  Well  known 
because  of  its  wide  applicability,  this  method  fails  to 
describe  many  interesting  problems  because  diffusion  theory 
is  completely  local  in  the  sense  that  each  cell  is  influenced 
only  by  its  nearest  neighbors.  One  consequence  of  this 
local  nature  is  that  the  shape  of  the  spatial  neutron  flux 
distribution  is  strongly  influenced  by  the  boundary  conditions. 
However,  it  is  precisely  at  the  boundary  that  this  solution 
technique  is  most  inaccurate.  For  spherical  systems  with  a 
radius  of  a  few  mean  free  paths,  all  points  in  the  system 
are  "close”  to  the  boundary  and  thus  diffusion  theory  is 
not  applicable. 

Clearly,  a  compromise  is  needed  between  the  simplicity 
of  diffusion  theory  and  the  accuracy  of  transport  theory. 


An  efficient  solution  technique  that  can  be  used  on  small  or 
medium  sized  computers  would  be  of  value  in  parametric 
studies  of  small,  fast,  supercritical  systems. 

Purpose  and  Criteria 

The  purpose  of  the  research  reported  here  was  to 
develop  a  technique  that  can  be  used  to  efficiently  calcu¬ 
late  the  neutron  flux  and  the  neutron  growth  rate  in  small 
spherically  symmetric,  multiplying  systems.  The  technique 
is  to  be  a  numerical  one  that  maximizes  the  use  of  analytical 
tools  and  approximations  in  order  to  minimize  the  numerical 
effort  and  expense.  Extreme  accuracy  is  sacrificed  in 
favor  of  efficiency,  but  the  strength  and  significance  of 
the  technique  will  lie  in  its  applicability  to  nonlocal 
problems  7-  problems  which  diffusion  theory  cannot  describe. 

Overview 

This  report  contains  five  sections  and  two  appendices. 

In  Section  II,  the  theoretical  background  of  some  semi- 
analytical  methods  that  influenced  the  development  of 
integral  equation  synthesis  (IES)  is  discussed.  The  extension 
of  these  methods  and  the  derivation  of  integral  equation 
synthesis  is  presented  in  Section  III.  In  Section  IV, 
benchmark  calculations  are  described  and  the  results  are 
analyzed.  Finally,  the  conclusions  and  recommendations  of 
this  study  are  presented  in  Section  V.  In  Appendix  A,  the 
numerical  methods  used  in  the  IES  calculations  are  described. 
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A  listing  of  pertinent  data  used  in  benchmark  calcula¬ 
tions  is  presented  in  Appendix  B. 
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II .  Background 


Early  attempts  to  solve  the  neutron  transport  equation 
were  based  on  approximations  that  allowed  analytical  solutions. 
As  computing  machinery  evolved  to  larger  and  faster  models, 
numerical  methods  also  evolved  to  brute  force  solutions  with 
fewer  and  fewer  approximations.  However,  many  of  the  early 
techniques  were  surprisingly  accurate  and  can  be  easily 
extended  to  describe  very  complex  problems  today. 

The  discussion  of  methods  and  techniques  of  neutron 
transport  included  in  this  section  is  certainly  not  intended 
to  be  a  comprehensive  survey  of  such  methods.  Rather,  the 
purpose  here  is  to  briefly  review  only  those  techniques  that 
provided  a  foundation  for  and  influenced  the  development  of 
integral  equation  flux  synthesis.  Some  of  the  basic  ideas 
of  synthesis  methods  were  detailed  over  three  decades  ago 
by  Serber  and  Wilson. 

Serber-Wilson  Method 

The  Serber-Wilson  method  (Ref  19)  is  an  approximate 
technique  for  finding  critical  masses  and  multiplication 
rates  for  two-media,  spherical  systems.  Developed  inde¬ 
pendently  by  Serber  and  Wilson  in  1945,  this  method 
approximates  the  neutron  flux  in  each  medium  with  asymp¬ 
totic  diffusion  theory  solutions.  The  integral  form  of 
the  neutron  transport  equation  is  required  to  be  satisfied 
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at  the  center  of  the  system.  This  requirement  provides 
the  constraint  necessary  to  define  a  matching  coefficient 
for  the  two  pieces  of  diffusion  theory  solutions. 

If  the  system  is  divided  into  a  core  region  and  a 
tamper  or  reflector  region,  then  diffusion  theory  can  be 
used  to  approximate  the  neturon  distribution  in  each  region. 
In  the  core  region  or  active  region,  the  neutron  density 
is  taken  to  be  of  sinusoidal  form;  in  the  reflector  region, 
an  exponential  is  used  for  the  neutron  density.  The  true 
distribution  shows  a  transition  effect  near  the  core¬ 
reflector  interface  where  the  density  drops  rapidly  in 
passing  from  the  core  to  the  reflector.  Serber  represented 
this  by  allowing  a  discontinuity  in  neutron  density  at  the 
interface.  The  magnitude  of  the  discontinuity  can  be 
determined  from  conservation  of  neutrons;  in  a  critical 
assembly  the  rate  of  neutron  production  in  the  core  must 
exactly  equal  the  rate  of  neutron  absorption  in  the  reflector 
(assumed  infinitely  thick).  Once  the  neutron  density  is 
fixed  in  this  way,  the  critical  radius  can  be  determined 
by  requiring  that  the  integral  equation  which  governs  the 
diffusion  of  neutrons  be  satisfied  at  r  ■  0 

The  striking  feature  of  the  Serber-Wilson  method  is 
that,  in  the  resulting  equation  to  be  solved  for  the  critical 
radius,  R  ,  one  side  of  the  equation  contains  only  core 
constants  and  the  other  side  of  the  equation  only  reflector 
constants.  This  obviously  leads  to  a  very  simple  graphical 
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solution.  The  extension  of  the  method  to  non-critical 
assemblies  is  straightforward  since  a  multiplying  system 
is  equivalent  to  a  critical  one  with  time  absorption 
(Ref  8:78). 

Deficiencies  in  the  Serber-Wilson  approach  are  three¬ 
fold:  first,  the  method  cannot  handle  systems  so  small  that 

the  diffusion  approximation  introduces  unacceptable  error; 
second,  although  the  method  can  be  extended  to  handle  three 
material  systems,  it  is  difficult  to  treat  multi-material 
systems  with  varying  densities;  third,  the  method  cannot 
easily  describe  the  effect  of  the  distribution  of  neutron 
energies  in  finite  systems  (Ref  6:275). 

In  Feynman's  method,  which  is  described  next,  special 
attention  is  paid  to  those  problems  arising  from  the  fact 
that  neutrons  of  different  velocities  have  different 
properties . 

Feynman's  Method 

Like  the  Serber-Wilson  method,  Feynman's  method 
(Ref  8)  is  an  approximate  technique  for  calculating  critical 
sizes  and  multiplication  rates  of  spherical,  active  cores 
surrounded  by  reflectors.  The  basis  of  Feynman’s  method 
is  the  integral  equation  for  A(x,v)  ,  the  density  of 
absorptions  of  neutrons  at  position  5t  in  the  core,  at 
speed  v  ,  per  unit  range  of  v  ,  Feynman  also  defines  the 
term  S(v'-*-v)dv  as  the  number  of  neutrons  (isotropic) 
that  emerge  in  velcoity  range  dv  when  a  neutron  of  velocity 
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v"  makes  a  collision  in  the  core.  He  further  defines  a 


kernel  P(x'->x,v)  as  the  density  of  absorptions  at  x  when 
one  neutron  is  released  isotropically  at  x'  with  velocity 
v  .  The  variable  v  merely  specifies  the  constants  to  be 
used  in  calculating  this  kernel.  Then  A(x,v)  satisfies 
the  following  integral  equation 

A(x,v)  =  Jdx'Jdv'  P(x+x,v)S(v'->-v)A(x%v')  (2.1) 


Now  a  form  for  the  kernel  P  can  be  obtained  by 
considering  a  simple  one-velocity  problem.  The  scalar  flux 
eigenfunctions  at  velocity  v  ,  <|>n(x,v)  satisfy  a  one- 

velocity  integral  equation 


vn (x,v) 


(2.2) 


where  kR(v)  are  the  associated  eigenvalues.  Feynman 
then  expands  the  kernel  of  Eqn  (2.2)  as  a  bilinear  series 
of  its  eigenfunctions. 


P(x'-«c,v) 


l  k»w 


*n (x,v)ipn(x"  ,v) 
Jdx  i|>£(x,v) 


(2.3) 


The  crux  of  Feynman’s  method  is  to  approximate  the  kernel  P 
in  various  ways  that  allow  a  simple  solution  to  Eqn  (2.1). 
The  first  lower  approximation  (so  named  because  it  under¬ 
estimates  the  critical  radius)  is  to  replace  all  the  kn(v) 
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in  Eqn  (2.3)  by  kQ(v)  so  that  then 
P(x'-*x,v)  =  kQ(v)  S(x'-x) 

and  after  substitution,  Eqn  (2.1)  simplifies  to 


A(x,v)  =  kQ  (v)  Jdv'  S(v'->v)  A(x,v') 


Similarly  the  first  upper  approximation  (it  overestimates 
the  critical  radius)  is  to  set  all  kn  (v)  =  0  except 
kQ(v)  .  Then  Eqn  (2.3)  becomes 


P(x'-»x,v) 


kQ(v) 


<(»0(x%v)  K»0(x,v) 
/dx  ^(x,v) 


Substituting  Eqn  (2.6)  into  Eqn  (2.1)  gives 


A(x,v) 


kQ(v) 


^0(x,v) 
/dx  iPq(x,v) 


An  implicit  assumption  of  Feynman's  approach  is  that 
the  one-group  problem  has  been  solved  for  every  velocity 
group.  In  practice,  he  used  the  diffusion  theory  solution 
for  the  i|»n(x,v)  ,  i.e.,  the  eigenfunction  <J>0  was  always 

approximated  by  -s - (Ref  9:163).  The  method  itself 
does  not  require  any  particular  approach  for  solving  the 
one-group  problem.  However,  Feynman's  method  does  require 


(2.4) 


(2.5) 


(2.6) 


(2.7) 
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that  the  spherical  system  under  consideration  be  divided 
into  a  core  and  a  reflector,  each  with  homogeneous  composi¬ 
tion.  If  this  is  not  the  case,  the  method  breaks  down.  If 
the  composition  is  nearly  uniform,  a  clumsy  perturbation 
treatment  can  still  be  made,  but  the  main  advantages  of  the 
method  are  lost  (Ref  8). 

Clearly.  Feynman's  method  gives  insight  into  an 
efficient  technique  for  determining  neutron  flux.  Rather 
than  solving  a  two-dimensional  problem  (one  spatial  dimension 
and  one  energy  dimension),  Feynman's  approximation  allows 
one  to  solve  two  one-dimensional  problems  successively. 
Although  this  is  strictly  possible  only  when  the  kernel 
is  actually  separable,  the  idea  of  Feynman's  method  is 
analogous  to  modern  flux  synthesis  methods. 


Flux  Synthesis  Methods 

Modern  flux  synthesis  methods  are  attempts  to  solve 
multidimensional  problems  by  first  solving  a  sequence  of 
simpler,  usually  one  dimensional  problems.  The  reconstructed 
multidimensional  solution,  although  an  approximation,  often 
contains  detail  which  is  impossible  or  economically  imprac¬ 
tical  to  achieve  with  a  direct  solution  (Ref  20). 

Early  flux  synthesis  methods  were  directed  toward 
solving  complex  reactor  problems  and  invariably  used 
diffusion  equations  and  thermal  energy  spectra.  Lancefield 
(Ref  13)  extended  space-energy  flux  synthesis  methods  by 
avoiding  the  diffusion  approximation  and  retaining  the 


governing  transport  equation.  He  also  introduced  several 
refinements  to  the  basic  space-energy  flux  synthesis  method. 
These  included:  leaving  both  the  space/angle  and  energy 
dependence  of  the  trial  function  to  be  determined  by  the 
variational  principle,  incorporating  discontinuous  trial 
functions,  and  the  use  of  a  new  variational  principle  for 
criticality  problems  that  leads  to  estimates  of  homogeneous 
functionals  of  the  unknown  flux. 

Lancefield’s  approach,  however,  is  still  plagued  with 
the  problem  of  choosing  appropriate  trial  energy  functions, 
<|>n(E)  .  One  customarily  chooses  the  <}>n(E)  to  represent 

infinite-media  spectra  characteristic  of  the  sub-regions  of 
the  system.  But  this  choice  is  restricted  to  large  systems 
where  such  solutions  dominate  in  some  regions.  For  small, 
highly  enriched  fuel  systems,  it  seems  likely  that  such  a 
choice  would  lead  to  considerable  errors  in  both  the  computed 
reaction  rates  and  the  calculated  energy  spectra  (Ref  21). 
Indeed,  Lancefield’s  results  for  a  two  region  (core  and 
reflector)  fast  reactor  gave  good  results  for  multiplication 
factor,  k  ,  but  the  errors  in  the  reflector  flux  were 
intolerably  large.  To  reduce  these  errors,  Lancefield  used 
different  spectra  and  weight  functions  in  different  regions. 
However,  the  spectra  and  weight  functions  were  obtained  from 
the  known  solution.  To  avoid  reliance  on  the  known  solution, 
an  iterative  scheme  was  devised,  but  Lancefield  himself 
admits  the  whole  procedure  became  so  complicated  and  lengthy 
that  any  practical  applicability  became  questionable. 


11 


Nearly  simultaneous  with  Lancefield's  work,  Neuhold 
and  Ott  (Ref  16)  improved  the  space-energy  synthesis  approach 
by  employing  reaction  rate  weighting,  by  using  realistic 
trial  functions,  and  by  deriving  a  more  general  analytical 
solution  for  the  synthesis  equations  which  includes  the 
use  of  complex  buckling.  These  improvements  led  to  error 
reductions  (compared  to  previous  space-energy  synthesis 
versions)  of  a  factor  of  100  in  multiplication  rate,  k  , 
and  a  factor  of  20  in  integral  quantities  sensitive  to  the 
non-separability  of  space  and  energy.  Unfortunately,  the 
applicability  of  these  improvements  is  limited  to  systems 
in  which  diffusion  theory  provides  an  adequate  model. 

Cockayne  (Ref  4)  extended  the  work  of  Neuhold  and  Ott 
by  examining  the  choice  of  trial  spectra.  He  found  that 
"realized"  spectra,  i.e.,  spectra  actually  existing  at  some 
point  in  the  system,  generally  gave  more  accurate  results 
than  a  spectrum  averaged  over  a  region.  Cockayne  also 
developed  a  procedure  that  allows  calculation  of  a  "realized" 
spectrum  in  any  transition  region. 

In  summary,  flux  synthesis  methods  offer  considerable 
potential  relative  to  reducing  computational  effort  for 
reactor  analysis.  However,  the  various  approaches  suffer 
from  lack  of  wide  applicability.  Most  of  these  approaches 
are  based  on  diffusion  theory  and  thus  the  applicability 
is  restricted  to  problems  in  which  the  system  is  large 
compared  to  the  operative  mean  free  paths.  Other  forms  of 


flux  synthesis  require  accurate  knowledge  of  the  energy 
spectrum  which  can  only  be  obtained  through  complex  iterative 
schemes.  The  concept  of  space -energy  flux  synthesis  is 
fundamental  to  the  new  method  described  in  the  next  section. 
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III.  Development  of  Integral  Equation  Synthesis 


In  this  section  we  develop  the  equations  that  comprise 
the  integral  equation  synthesis  method.  Basically,  the 
goal  is  to  expand  the  two  dimensional  flux,  i.e.,  radius 
and  energy,  in  terms  of  products  of  one  dimensional,  sepa¬ 
rated,  trial  eigenfunctions.  The  one  dimensional  problems 
are  computationally  easy  and  economical  to  solve,  and  the 
accuracy  of  the  final  solution  can  be  adjusted  by  using 
more  trial  eigenfunctions. 

First,  the  integral  form  of  the  transport  equation  is 
specialized  and  simplified  for  a  one  dimensional,  mono- 
energetic,  time  independent  problem.  Analogous  to  this 
spatial  integral  equation  and  based  on  neutron  conservation, 
we  next  develop  the  energy  integral  equation.  Finally, 
the  synthesis  procedure  is  described  whereby  the  expansion 
coefficients  are  determined  and  the  final  solution  is 
constructed . 

Spatial  Integral  Equation 

The  basis  of  the  integral  equation  synthesis  method 
is  the  integral  form  of  the  neutron  transport  equation 
(Ref  6:27) 

,P 

-/  E(E,p'jdp ' 

*(r,fl,E,t)  =  JdV'  - - — -  JdE'  $(r ',$',E',t') 

4.  1* 

E(r',E')C(r',E')fr-(E'-E;  n'-fl)  C3>1) 
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where 


$(r,ft,E,t)  =  neutron  flux  at  point  r  , 

energy  E  ,  solid  angle  ft  , 
at  time  t 


E(r',E') 


P 

C(r-,EO 


total  macroscopic  cross  section 
at  (r',E') 


average  number  of  secondary  neutrons 
emitted  per  collision  at  (r^E") 

probability  that  a  neutron  entering 
a  collision  with  energy  E'  and 
direction  ft'  will  emerge  with 
energy  E  and  direction  ft 


Several  simplifications  and  assumptions  are  now  made 
in  order  to  reduce  the  number  of  independent  variables  in 
the  flux  that  is  to  be  computed.  First,  the  angular  varia¬ 
tion  of  the  flux  is  eliminated  by  assuming  isotropic  scatter 
and  integrating  both  sides  of  Eqn  C3.1)  over  all  solid  angle. 
From  now  on,  the  all  angle  or  scalar  flux  is  being  computed. 
It  will  become  apparent,  in  the  development  below  that  this 
simplification  is  not  essential  to  the  method,  but  is  a 
significant  convenience. 
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Secondly,  the  time  dependent  problem  can  be  transformed 
into  a  stationary  problem  (Ref  8:78).  We  assume  the  time 
dependence  of  the  flux  in  a  multiplying  system  is  given  by 
an  exponential  growth  rate  a  .  Thus, 

*(r,E,t)  =  *(r,E)  eat  (3.2) 

and 

*(r',E',t  -  )  =  t(r',E')  eat  e  (3.3) 

Substitution  of  Eqns  (3.2)  and  (3.3)  into  (3.1)  transforms 
the  integral  equation  into  a  time  independent  problem. 

A  further  simplification  is  made  by  restricting 
applications  to  one-dimensional,  spherically  symmetric 
systems.  The  specialized  integral  equation  is  now 

P 

-J  £(p',E)dp' 

o 

*(r,E)  =  JdV'  - -  fdE'  *(r ' ,E") 

4  TT  p  2 

- 

e  V^EJ  E(r',E')C(r',E')fr-(E^E)  (3.4) 
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Note  that  the  time  retardation  term,  exp  [-ap/v(E)]  ,  can 

be  incorporated  into  the  exponential  attenuation  term  by 
defining  a  new  cross  section 

Z*(p,F)  =  l(p,E)  +  (3.5) 

We  now  write  the  integral  equation  for  one  energy  group. 

-JPE*(p')dp' 

o 

x«Kr)  =  JdV'  - -  <Kr')E(r')C(r')  (3.6) 

4  7T  p  ^ 

where 

\|i(r)  =  one  group  spatial  flux  distribution 

i*(p')  =  i(p')  +  -2- 

A  =  eigenvalue 

Equation  (3.6)  represents  an  auxiliary  problem  which  is  an 
eigenvalue  problem.  It  can  be  solved  numerically,  such  as 
by  the  power  method  (Ref  7:291),  for  the  spatial  flux 
distribution  i|>(r)  .  In  addition,  the  lowest  order  eigen¬ 
value  of  Eqn  (3.6)  can  be  made  equal  to  unity  (i.e.,  an 
exactly  critical  system)  by  iterative  adjustments  to  I* 
as  per  Eqn  (3.5).  When  the  eigenvalue  is  made  equal  to  one, 

Eqn  (3.5)  can  be  solved  for  the  one  group  multiplication 
rate,  o 
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However,  Eqn  (3.6)  has  an  infinite  number  of  solutions 
corresponding  to  higher  order  eigenvalues  and  associated 
eigenfunctions.  There  will  be  a  lowest  A  ,  called  A^  , 
for  which  the  eigenfunction  ^  is  everywhere  positive. 

The  eigenfunctions  corresponding  to  higher  order  values  of 
A  (x  being  associated  with  i|>n  )  will  have  one  or  more 
nodal  points.  In  other  words,  the  neutron  density  will  be 
positive  in  some  regions  and  negative  in  others.  Feynman 
(Ref  8:8)  interprets  negative  neutron  densities  as  defi¬ 
ciencies  below  some  positive  neutron  density.  He  further 
states  that  physical  reasoning  can  be  used  in  interpreting 
the  integral  equation  if  negative  neutrons  are  thought  of 
as  actual  particles  whose  presence  in  a  region  can  cancel 
the  presence  of  an  equal  number  of  positive  neutrons. 
Clearly,  increasing  the  number  of  nodal  points  in  4«  will 
decrease  the  overall  significance  or  contribution  from 
that  i|>n  ,  This  is  so  because  leakage  of  neutrons  of 
opposite  sign  into  one  another  will  become  more  rapid  as 
the  eigenfunction  oscillates  more  rapidly. 

It  should  be  noted  here  that  the  eigenfunctions,  ipn  , 
are  not  mutually  orthogonal  since  the  kernel  of  Eqn  (3.6) 
is  not  symmetric.  Furthermore,  there  is  no  guarantee  that 
the  functions  form  a  complete  set,  even  though  this  is 

certainly  true  for  cases  of  physical  interest  (Ref  8). 

But  as  Kaplan  (Ref  12)  argues,  "completeness  is  academic  and 
orthogonality  is  only  a  minor  convenience.”  This  is 
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especially  true  for  the  case  of  an  approximate,  numerical 
solution.  Indeed,  at  most  only  a  few  of  the  set  of  trial 
eigenfunctions  will  be  used  in  the  final  solution.  The  only 
criterion,  then,  for  these  eigenfunctions  is  the  "goodness" 
of  the  final  approximation,  i.e.,  the  accuracy  of  the  final 
solution. 

We  now  turn  our  attention  to  the  energy  dependence  of 
the  neutron  transport  equation.  The  next  goal  is  to  develop 
an  energy  integral  equation  that  is  analogous  to  the  spatial 
integral  equation. 

Energy  Integral  Equation 

An  analogous  energy  integral  equation  can  also  be 
derived  from  Eqn  (3.4).  The  procedure  is  to  assume  spatially 
invariant  material  properties  and  a  system  of  infinite 
extent.  These  assumptions  allow  the  spatial  integrations 
of  Eqn  (3.4)  to  be  performed  analytically  resulting  in 
Eqn  (3.11).  However,  the  development  of  the  energy  integral 
equation  presented  below  is  based  on  conversation  of  neutrons 
and  serves  to  give  a  better  physical  interpretation  of  the 
energy  equation. 

Let  n(E)  be  the  energy  dependent  neutron  number 
density  in  units  of  neutrons  per  unit  volume  per  unit  energy. 
Then  the  time  rate  of  change  of  the  neutron  number  density 
is  just  equal  to  the  gains  minus  the  losses.  That  is. 


* 
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gains  -  losses 


(3.7) 


3n(E)  _ 

3t 

We  assume  a  system  of  infinite  extent  with  a  spatially 
constant  flux.  The  contribution  of  n(E)  ,  summed  for  each 
energy  bin  dE '  ,  is  given  by 


gains  =  JdE'  n(E')v(E')S(E>)C(E>)f(E'^E)  (3.8) 

where 


n(E  ) 
v(E') 
E(E') 


C  (E  ') 


f  (E'h-E) 


neutron  number  density  at  energy  E' 

neutron  speed  at  energy  E' 

macroscopic  total  cross  section  at 
energy  E' 

average  number  of  secondary  neutrons 
emitted  when  a  neutron  of  energy  E' 
undergoes  a  collision 

probability  that  a  neutron  entering  a 
collision  with  energy  E'  will  emerge 
with  energy  E 


Equation  (3.8)  is  nearly  self-explanatory.  The  number 
density  times  the  speed  is  just  neutron  flux.  The  product 
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of  flux  and  cross  section  is  the  reaction  rate.  For  each 
reaction  there  are  C(E')  neutrons  emitted  with  energy 
distribution  given  by  f(E'->E)  .  In  order  to  obtain  the 

entire  gain  to  n(E)  ,  we  simply  sum  or  integrate  over  all 
energies  E" 

The  losses  term  in  Eqn  (3.7)  is  even  more  straight¬ 
forward  than  the  gains.  Since  we  have  assumed  an  infinite 
system,  there  is  no  leakage.  The  only  loss  is  due  to  reactions 
or  collisions  that  occur  at  energy  E  .  This  is  simply  the 
product  of  flux  and  total  cross  section. 

losses  =  n(E)  v(E)  (E)  (3.9) 

Equation  (3.7)  can  now  be  written  as 

an^  =  fdE'  n(E')  v(E')  l(E')  C(E')  f(E'->E) 

3t 

-  n (E)  v (E)  (E)  (3.10) 

But  again,  as  in  Eqn  (3.2),  we  claim  to  know  the  time 
dependence  of  our  system.  In  this  case,  the  exponential 
growth  rate  is  given  by  aw  since  we  have  assumed  an 
infinite  system.  Taking  the  partial  derivative  in  Eqn  (3.10) 
and  rearranging,  we  have  the  energy  integral  equation. 
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[a  +  v (E) E (E) ]  n(E) 

oo 


jdE'  n(E')v(E')Z(E')C(E')f(E'-*-E)  (3.11) 

Note  that  Eqn  (3.11)  is  another  eignevalue  problem.  The 
lowest  order  solution  n^(E)  ,  corresponds  to  the  energy- 

distribution  in  a  homogeneous,  infinite  medium.  Just  as  in 
the  spatial  case,  there  are  an  infinite  number  of  solutions 
corresponding  to  higher  order  eigenvalues  and  associated 
eigenfunctions.  The  lowest  order  solution,  n^  (E)  ,  will 

be  everywhere  positive  and  higher  order  eigenfunctions  will 
have  one  or  more  nodal  points.  The  energy  eigenfunctions 
are  also  not  mutually  orthogonal,  although  the  adjoint  or 
left  eigenfunctions  are  orthogonal  to  the  right  eigen¬ 
functions.  Completeness  of  the  energy  set  of  functions 
is  again  academic  since  at  most  only  a  few  eigenfunctions 
will  be  used  to  approximate  the  final  solution. 

We  now  have  methods  for  determining  the  spatial  flux 
independent  of  energy,  and  also  for  determining  the  energy 
distribution  independent  of  spatial  location.  The  next 
step  is  to  combine  or  synthesize  the  solutions  to  the 
separated  problem.  The  final  solution,  a  judicious  combina¬ 
tion  of  the  separate  solutions,  will  be  the  integral 
equation  synthesis  approximation  to  the  two  dimensional, 
i.e.,  space  and  energy,  non-separable  problem. 
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Synthesis  of  Spatial  and  Energy  Solutions 

Assuming  the  spatial  and  energy  eigenfunctions  form  a 
complete  set,  the  neutron  flux  can  be  expanded  in  terms  of 
these  trial  eigenfunctions: 


’(r ,E)  =  l  |  a.,  (r)  n,  (E) 

i  i  1J  1  J 


(3 


Since  the  sets  of  spatial  and  energy  trial  functions 
are  known,  a  good  approximation  to  the  flux  should  be 
available  once  a  few  of  the  expansion  coeff icienct ,  a^j  , 
are  obtained. 

In  order  to  simplify  the  bookkeeping  notations,  we 
arrange  the  trial  solutions  in  an  arbitrary  but  consistent 
order  so  that  trial  solution  number  one  is  t|>^(r)  n^E)  » 
trial  solution  number  two  is  <|i2(r)  n^  (E)  ,  trial  solution 

number  three  is  i^(r)  n2(E)  ,  trial  solution  four  is 

<J»2  (r)  n2(E)  ,  etc. 

The  standard  methods  of  perturbation  theory  (Ref  14:413) 
are  used  to  obtain  the  expansion  coefficients.  That  is, 
we  first  operate  on  the  trial  solution  with  the  exact,  two- 
dimensional  kernel  of  Eqn  (3.4).  In  other  words,  replace 
$(r',E')  on  the  right  hand  side  of  Eqn  (3.4)  with  the 
known  i|^(r')  n1(E')  .  It  is  important  to  note  that 

Eqn  (3.4)  is  no  longer  an  eigenvalue  problem;  it  simply 
represents  an  integration  that  must  be  performed  to  obtain 
the  flux  *(r,E)  .  This  integration  can  be  performed 


12) 
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numerically  quite  quickly  since  no  iteration  or  simultaneous 
solutions  to  equations  are  required. 


Let  the  result  of  the  integration  be  (r^E)  ,  the 

subscript  denoting  that  the  first  trial  solution  was  used 
on  the  right  hand  side  of  Eqn  (3.4).  We  can  then  use  adjoint 
weighting  to  expand  (r^E)  in  terms  of  all  the  trial 

eigenfunctions  ^  n^  and  let  the  expansion  coefficients 
by  Aln  .  The  first  subscript  on  A  denotes  that  X^ 
is  being  expanded  and  the  second  subscript  identifies  the 
coefficient  as  belonging  to  the  n1*1  trial  solution. 

That  is , 


N 


max 


C1  (tjE)  I  Aln  ^  nj 

n=l  J 


1,2,,.. 


(3.13) 


and  clearly 


//X  (r1E)1|)|(r)n{(E)drdE 

Aln  =  - — - 1 - L -  (3.14) 

jf  nt  n^.  drdE 

where  the  dagger  terms  are  the  left  or  adjoint  eigen¬ 
functions  found  from  Eqn  (3.6)  for  spatial  functions  and 
from  Eqn  (3.11)  for  energy  functions.  The  denominator 
in  Eqn  (3.14)  can  easily  be  made  unity  with  proper  normali¬ 
zation  of  the  separate  eigenfunctions. 
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t 


Similarly,  the  integration  indicated  in  Eqn  (3.4)  is 
performed  after  substituting  ^2nl  *  ^ln2  ’  ^2n2  *  etc' 

for  $(r',E'j  on  the  right  hand  side  of  Eqn  (3.4).  The 


results  of  these  integrations  are  labelled  X2  ,  X^  , 

X^  ,  respectively,  and  these  in  turn  are  expanded  in  terms 
of  the  basis  set.  The  expansion  coefficients, 

defined  by  Eqn  (3.14),  form  a  transformation  matrix  [A  1 
This  matrix  transforms  from  the  i/^(r)  n^  (E)  basis  to  the 
X(x^E)  basis. 

Note  that,  if  the  trial  eigenfunctions  <|^  n^.  were 
exact  eigenfunctions  of  the  kernel  of  Eqn  (3.4),  i.e.,  of 
the  fully  two-dimensional,  non-separated  kernel,  then  the 
matrix  [A  ]  would  be  the  identity  matrix.  Thus,  the  mag¬ 
nitude  of  the  off-diagonal  elements  of  [Amn]  give  a  measure 
of  the  validity  of  the  separability  assumptions  used  to  obtain 
Eqns  (3.6)  and  (3.11).  In  addition,  if  the  choice  of  trial 
eigenfunctions  is  reasonably  accurate,  then  the  magnitude 
of  the  off-diagonal  elements  should  decrease  the  further  one 
goes  away  from  the  main  diagonal. 

The  next  step  in  the  synthesis  of  the  spatial  and 
energy  solutions  is  to  use  the  power  method  (Ref  7:291) 
to  obtain  the  largest  eigenvalue  and  associated  eigenvector 
of  the  transformation  matrix  [A  1  .  The  elements  of  the 

lowest  order  eigenvector  are  just  the  expansion  coefficients 
a^j  that  we  seek  for  Eqn  (3.12). 
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Now  that  a  new  and  much  more  accurate  approximation 
to  the  flux  is  available  through  Eqn  (3.12),  we  can  use 
this  flux  one  more  time  in  Eqn  (3.4)  to  perform  the  right 
hand  side  integration.  This  final  integration  will  yield 
a  more  accurate  eigenvalue  which  can  then  be  used  to  compute 
the  final  growth  rate,  a  .  We  then  have  solved  for  the 
flux  as  a  function  of  radius  and  of  energy, and  have  also 
determined  the  growth  rate  for  the  system. 

As  pointed  out  in  Section  I,  the  objective  of  this 
research  was  to  develop  an  efficient  numerical  technique 
for  computing  neutron  flux  and  growth  rate  in  specialized 
systems.  However,  the  entire  development  in  this  section 
has  been  in  terms  of  continuous  functions  and  continuous 
variables  which  are  obviously  not  suitable  for  direct 
transfer  to  digital  computers.  The  computer  encoding, 
numerical  techniques  and  approximations  are  certainly  not 
a  trivial  portion  of  any  solution  algorithm.  Indeed,  in 
the  case  of  integral  equation  synthesis,  a  large  part  of 
the  numerical  efficiency  can  be  attributed  to  a  new  and 
unique  method  of  volume  integration  on  a  sphere  (.Ref  17). 
However,  the  particular  numerical  methods  remain  secondary 
to  the  essential  ideas  of  integral  equation  synthesis. 

The  development  and  description  of  numerical  techniques  used 
in  integral  equation  synthesis  can  be  found  in  Appendix  A. 
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In  the  next  section,  integral  equation  synthesis  is 
used  to  obtain  neutron  flux  and  growth  rate  for  small, 
spherical  systems.  These  results  are  compared  to  those 
obtainable  with  more  conventional  (and  expensive) 
calculational  schemes. 


IV.  Results  and  Discussion 

The  integral  equation  synthesis  (IES)  approach  was 
converted  into  a  Fortran  computer  program  by  the  author 
and  used  to  compute  the  neutron  flux  and  growth  rate  in 
some  simple  systems.  In  this  section  the  sample  problems 
are  described,  as  well  as  the  results  obtainable  with  this 
method . 

Jezebel  Calculations 

Jezebel  is  a  bare  plutonium  sphere,  measured  to  be 
just  critical  with  a  mass  of  16.6  kg  at  a  density  of 

15. 8  g/cm3,  corresponding  to  a  radius  of  6.4  cm.  The 
isotropic  composition  of  the  plutonium  is  94.1  wt  %  Pu-239, 

4.8  wt  %  Pu-240,  and  1.0  wt  %  gallium  (Ref  5). 

The  benchmark  for  IES  calculations  is  the  experimental 
data  reported  in  LA-3529  (Ref  5).  In  addition,  flux  and 
energy  distributions  were  obtained  with  the  1977  version 
of  the  Los  Alamos  DTK  code  using  their  standard  13  energy 
group  cross  section  set.  The  results  from  the  DTK  code 
differ  insignificantly  from  those  quoted  in  LA-3529,  and 
so  these  sources  are  used  interchangeably  as  a  benchmark 
and  are  herein  referred  to  as  LASL  results. 

Cross  Sections  Used  in  Calculations 

The  cross  sections  used  in  the  IES  calculations  were 
a  16  group  set  collapsed  from  the  175  group  set  reported 


in  UCRL-50400  Vol .  16  (Ref  18).  The  collapsing  was  done  by 
Dr.  Walt  Webster  of  Lawrence  Livermore  Laboratory  (LLL) 
who  used  a  group  structure  and  weighting  spectrum  commonly 
used  at  LLL  for  small,  fast,  critical  assemblies.  The 
group  structure  and  associated  group  velocities  are  described 
in  Appendix  B. 

It  is  important  to  note  here  that  the  cross  sections 
used  in  the  IES  calculation  are  not  identical  to  those 
used  in  the  benchmark  calculation.  This  selection  was 
made  for  two  reasons.  First,  the  LLL  cross  sections  were 
readily  available  and  in  a  form  easily  incorporated  into  the 
IES  methodology.  Secondly,  and  more  importantly,  it  was 
decided  early  in  this  project  that  the  level  of  accuracy 
desired  with  the  IES  method  was  inconsistent  with  highly 
refined,  normalized  cross  section  sets  that  give  good 
results  in  only  one  specific  application.  Indeed,  the 
authors  of  the  benchmark  calculation  point  out  that  cross 
section  users  often  find  that  normalization  is  essential 
to  achieve  the  accuracy  required  of  their  calculations 
(Ref  11).  The  results  reported  for  the  IES  calculations 
are  typical  of  those  obtainable  with  any  reasonable  cross 
section  set. 

Numerical  Details 

To  perform  the  integrations  needed  in  the  integral 
equation  synthesis  method,  Gaussian  quadratures  are  used 
for  the  angular  integrations.  The  Gaussian  weights  and 
abscissas  are  given  in  Appendix  B. 
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The  spatial  mesh  used  in  the  IES  Jezebel  calculation 


consists  of  a  center  ball  and  17  shells  of  increasing 
thickness.  For  ease  of  calculation  and  numerical  integration, 
a  self -similar  mesh  is  chosen  such  that  the  outside  radius 
of  any  shell  is  given  by 


r 


i 


e 


r 


i-1 


(4.1) 


where 


r^  =  outside  radius  of  an  arbitrary  shell 

ri  ;L  =  outside  radius  of  next  inner  shell 

e  =  1.1571 

Growth  Rate  Results 

Integral  equation  synthesis  computed  the  Jezebel 
alpha  to  be  +0.39  generations  per  microsecond,  compared 
to  the  benchmark  value  of  -  0.65  generations  per  micro¬ 
second.  A  better  appreciation  of  the  accuracy  of  the  IES 
result  can  be  obtained  by  converting  the  alpha  values 
(which  should  be  exactly  zero  for  a  critical  assembly)  to 
multiplication  factor,  k  (which  is  exactly  one  for  a 
critical  assembly).  The  LASL  value  corresponds  to  a  k 
of  0.99992  whereas  the  IES  alpha  corresponds  to  a  k 
of  1.00005.  This  deviation  of  less  than  0.01%  is 
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indeed  negligible  in  view  of  uncertainties  associated  with 
cross  sections  and  other  material  properties  such  as 
isotopic  composition  and  density. 

Computed  Flux  Results 

The  spatial  flux  distributions  computed  by  the  IES 
method  and  by  the  LASL  code  are  shown  in  Fig.  2.  The  two 
curves  are  normalized  to  each  other  at  0.71  cm  since  the 
different  meshes  in  each  calculation  had  a  common  value  at 
this  radius.  This  normalization  is  convenient  and  yet, 
for  practical  purposes,  it  is  equivalent  to  normalization 
at  the  center  of  the  assembly;  i.e.,  at  zero  radius. 
Normalization  at  the  center  is  preferred  to  an  integral 
normalization  because,  in  a  highly  supercritical  assembly, 
the  central  region  drives  the  entire  system  so  that  errors 
near  the  center  are  much  more  significant  than  differences 
near  the  outer  edge. 

As  is  typical  of  IES  results,  the  method  underestimates 
the  magnitude  of  the  flux  in  the  outer  regions.  Note  that 
the  maximum  discrepancy  in  the  flux  at  the  outer  edg j  of 
the  assembly  is  only  about  8%.  This  is  most  likely  due  to 
an  approximation  used  in  performing  the  volume  integrations 
required  by  the  integral  equation  approach.  This  approxi¬ 
mation  is  not  essential  to  the  method,  but  does  significantly 
reduce  the  calculational  effort.  See  Appendix  A  for  a 
detailed  description  of  the  geodesic  coefficients  used  in 
volume  integrations. 
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Figure  2.  Spatial  Flux  Distribution  for  Jezebel  Assembly 
as  Calculated  by  the  LASL  Code  and  IES  Method 


The  IES  results  shown  in  Fig.  2  were  obtained  with 
only  two  trial  eigenfunctions.  That  is,  the  final  spatial 
flux  distribution  is  a  combination  of  Cr )  and  ^(r) 

only.  Plotted  in  Fig.  3  are  the  LASL  result,  the  IES 
result,  and  the  ^  only  result.  Note  that,  in  the 
Jezebel  case,  the  addition  of  another  trial  eigenfunction 
(IES  curve  compared  to  curve)  improves  the  result,  but 

not  dramatically.  Just  as  one  would  expect,  the  lowest 
order  trial  eigenfunction  carries  most  of  the  information 
and  additional  trial  eigenfunctions  represent  fine  tuning 
of  the  basically  accurate  first  approximation.  In  many 
practical  applications,  it  may  be  determined  that  the 
improvement  in  accuracy  obtained  by  introducing  more  trial 
eigenfunctions  is  not  necessary  or  not  worth  the  computa¬ 
tional  effort.  The  integral  equation  synthesis  method  is 
inherently  flexible  and  easily  adapted  to  various  accuracy 
requirements . 

The  computed  energy  distribution  of  the  Jezebel 
flux  at  the  outer  surface  is  shown  in  Fig.  4.  Also  plotted 
in  Fig.  4  is  the  LASL  leakage  spectrum  for  the  Jezebel 
device.  The  two  histograms  are  normalized  to  the  same 
value  at  the  peak  energy  groups  respectively  which  occur 
near  .5  MeV.  Although  the  group  structures  are  different, 
not  the  general  agreement  between  the  IES  results  and  the 
LASL  results.  Just  as  in  the  spatial  case,  the  final  IES 
energy  distribution  was  obtained  with  only  two  trial 
eigenfunctions,  n^(E)  and  ^(E)  .  But  contrary  to  the 


Effect  of  Additional  Trial  Eigenfunction 
on  the  Jezebel  Spatial  Flux  Distribution 


spatial  case,  the  addition  of  the  second  energy  eigenfunction 
made  only  negligible  changes  in  the  computed  spectrum. 
Obviously,  the  infinite  medium  spectrum  is  very  close  to 
the  Jezebel  spectrum.  In  other  words,  essentially  all  the 
spectral  information  is  carried  in  the  lowest  order  trial 
eigenfunction.  We  would  not  expect  such  a  fortuitous 
computation  of  energy  eigenfunctions  in  systems  of  variable 
density  or  non-homogeneous  composition.  In  any  event, 
when  additional  energy  eigenfunctions  do  not  significantly 
alter  the  final  spectrum,  it  is  clearly  uneconomical 
and  unnecessary  to  include  more  than  one  trial  function. 

Shown  below  in  Table  I  is  the  transformation  matrix 
computed  for  the  Jezebel  system.  Note  that  the  terms  of 
maximum  magnitude  occur  on  the  main  diagonal. 


Table  I 

Transformation  Matrix  from  Jezebel  Calculation 

Vi 

*2nl 

V2 

*2n2 

X1 

1.0358 

-0.0062 

0.0019 

-0.0004 

X2 

-0.0035 

0.5728 

-0.0003 

-0.0015 

X3 

0.1891 

-0.0011 

-0.6621 

-0.0218 

X4 

-0.0006 

0.1059 

-0.0156 

-0.5149 

i 


In  Table  II  below  are  the  expansion  coefficients 
computed  from  the  transformation  matrix  shown  in  Table  I. 
Clearly,  the  maximum  information  is  "contained"  in  the 


lowest  approximation;  i.e.,  in  •  It  is  also  clear 

from  the  relative  magnitude  of  the  coefficients  shown  in 
Table  II  that  n2  carries  little  information.  In 
comparison  the  second  spatial  trial  function,  ^  , 

contributes  about  10%  to  the  final  solution. 


Table 

Mixing  Coefficients  for 

II 

Jezebel  Flux  Synthesis 

Function 

Coefficient 

*lnl 

0.9938 

*2nl 

0.1107 

hn2 

-0.0076 

*2n2 

-0.0020 

Computation  Times 

A  Fortran  program  using  the  IES  method  has  been 
executed  on  an  IBM  360/75  digital  computer.  As  with  all 
calculations  of  this  type,  central  processing  unit  (CPU) 
times  are  very  problem  dependent.  The  Jezebel  calculation 


with  17  spatial  shells,  16  energy  groups,  and  4  trial 
eigenfunctions  (2  spatial  and  2  energy)  requires  about 
90  seconds  CPU  time.  The  same  problem  with  17  spatial 
shells,  4  energy  groups,  and  only  the  lowest  order  spatial 
and  energy  eigenfunctions  requires  about  10  seconds  CPU 
time. 

In  comparison,  ANISN  (an  anisotropic,  discrete 
ordinates  transport  code  designed  for  reactor  analysis) 
requires  about  five  minutes  of  CPU  time  on  a  CDC  7600 
computer  (Ref  22).  Other  very  sophisticated  codes  used 
at  LLL  to  treat  small  systems  require  about  ten  minutes 
of  CDC  7600  CPU  time  to  compute  the  Jezebel  flux  and 
growth  rate  (Ref  22). 

Supercritical  and  Non-Homogeneous  Systems 

Integral  equation  synthesis  has  been  used  to  calculate 
various  supercritical  systems  including  single-material, 
variable-density  assemblies  and  also  several -material , 
variable-density  systems.  Unfortunately,  it  is  difficult 
if  not  impossible  to  find  benchmark  calculations  of  this 
type  in  the  open  literature.  Yet  it  is  just  these  more 
sophisticated  systems  that  integral  equation  synthesis  is 
designed  to  calculate. 

A  validation  of  a  portion  of  the  IES  method  has  been 
made  by  solving  a  one  group,  supercritical  system  by 
several  methods.  A  hypothetical  assembly  consisting  of  a 
5  cm,  double  density,  bare,  plutonium  sphere  was  analyzed 


using  diffusion  theory,  a  P-1  model,  discrete  ordinates, 
the  end-point  method,  and  the  spatial  part  of  the  integral 
equation  synthesis  method.  Table  III  shows  the  computed 
values  for  the  growth  rate  a  ,  and  the  corresponding 
multiplication  factor,  k  (Ref  3). 


Table  III 

Comparison  of  Calculated  Alpha  Values 


Model 

Alpha 

(Sec" 1 ) 

k_ 

Diffusion 

2.424 

X 

108 

1.67 

P-1 

2.727 

X 

108 

1.82 

Discrete  Ordinates 

2.973 

X 

108 

1.96 

End-Point 

2.965 

X 

108 

1.96 

Integral  Equation  Synthesis 

2.987 

X 

108 

1.97 

For  this  problem  the  end-point  method  can  be  considered 
as  nearly  exact  (Ref  2:96)  and  hence  serves  as  the  bench¬ 
mark.  Notice  that  diffusion  theory  predicts  an  a  that 
differs  by  181  from  the  end-point  method.  Integral  equation 
synthesis,  on  the  other  hand,  predicts  an  a  within  0.07% 
of  the  end-point  value. 

Integral  equation  synthesis  has  also  been  applied  to 
multi-material,  variable  density,  supercritical  systems. 
Benchmark  results  are  not  available  for  these  assemblies; 


t- 
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however 

results 


it  is  anticipated  that  the  accuracy  of  the  IES 
would  be  consistent  with  the  accuracies  indicated 


above. 


V.  CONCLUSIONS  AND  RECOMMENDATIONS 


Conclusions 

Integral  equation  synthesis  has  been  shown  to  be  an 
accurate  and  efficient  technique  for  calculating  neutron 
flux  and  growth  rate  in  spherical  systems.  Because  an 
integral  approach  is  used,  the  numerical  technique  does 
not  suffer  from  convergence  problems  so  common  with  finite 
difference  schemes.  Because  the  trial  functions  are  com¬ 
puted  for  each  problem,  the  method  does  not  suffer  from 
the  occasional  anomalous  results  of  usual  synthesis  methods 
(Ref  15)  where  the  trial  functions  are  chosen  a  priori. 

And,  because  the  method  uses  an  integral  approach,  the 
boundary  conditions  are  automatically  and  rigorously 
satisfied  (Ref  10). 

The  computational  efficiency  of  the  integral  synthesis 
method  arises  not  only  from  the  method  itself,  but  also 
from  the  coarser  mesh  that  can  be  used  in  most  applications. 
The  use  of  a  self-similar  mesh  not  only  simplifies  the 
calculations,  but  also  serves  to  give  maximum  spatial 
resolution  of  the  flux  near  the  center  of  the  system.  It 
is  precisely  this  region  that  drives  the  entire  assembly 
so  that  the  fine  resolution  is  well-placed.  On  the  other 
hand,  the  regions  near  the  outside  edge  of  a  highly  super¬ 
critical  assembly  contribute  little  to  the  center  flux 
and  overall  growth  rate.  The  self-similar  mesh  is 


42 


justifiably  very  coarse  near  the  outer  edge  of  the  system. 

The  mixing  coefficients  or  the  actual  synthesis  of 
the  separate  space  and  energy  solutions  has  been  shown  to 
be  unnecessary  for  acceptable  accuracy  in  Jezebel  calcula¬ 
tions.  Clearly,  the  synthesis  procedure  is  not  necessary 
in  problems  that  are  nearly  separable  in  terms  of  space  and 
energy.  Fortunately,  the  method  itself  gives  an  indication 
of  the  accuracy  of  the  separability  assumption.  The  trans¬ 
formation  matrix,  whose  lowest  order  eigenvector  produces 
the  mixing  coefficients,  provides  a  measure  of  the  separa¬ 
bility  through  the  magnitude  of  the  off-diagonal  elements. 

If  these  elements  are  large,  say  within  an  order  of  magnitude 
of  the  diagonal  elements,  then  the  synthesis  procedure  is 
probably  required  because  of  non-separability.  However,  if 
the  synthesis  procedure  is  not  necessary,  the  integral 
method  is  an  even  more  efficient  technique  for  computing 
neutron  flux  in  small,  supercritical  systems.  A  unique 
feature  of  integral  equation  synthesis  is  that  an  interim 
calculational  produce  gives  an  indication  of  the  accuracy 
of  the  final  product. 

The  wide  applicability  of  this  transport  technique 
cannot  be  overemphasized.  The  method  is  only  slightly 
more  complicated  than  diffusion  theory,  yet  it  handles  a 
much  broader  range  of  systems.  There  are  no  "size  compared 
to  mean  free  path"  restrictions  as  with  diffusion  theory. 

The  method  is  applicable  to  multi -material  systems  with 


varying  densities  and/or  systems  with  center  voids.  Yet 
the  method  accounts  for  the  nonlocal  nature  of  neutron 
interactions  in  a  multiplying  system.  The  applicability 
of  integral  equation  synthesis  is  further  enhanced  because 
of  its  modest  computing  cost  in  terms  of  storage  and 
CPU  time. 

Recommendations 

Two  extensions  of  this  research  and  integral  equation 
synthesis  are  recommended.  First,  the  method  should  be 
extended  to  allow  anisotropic  scatter  and  to  obtain  the 
angular  distribution  of  the  neutrons.  As  it  is  well  known, 
the  angular  distribution  of  the  neutrons  can  be  reproduced 
from  the  knowledge  of  the  flux  by  a  simple  quadrature 
(Ref  10).  This  technique  is  especially  straightforward 
the  the  numerical  form  of  the  IES  method  since  Gaussian 
quadrature  is  already  used  to  perform  the  angular  portion 
of  the  volume  integrations.  Instead  of  summing  the  contri¬ 
butions  of  all  angular  rays,  one  has  to  simply  retain  and 
store  the  contribution  of  each  angular  ray.  In  other 
words,  the  angular  dependence  of  the  neutron  flux  is 
already  being  computed  and  all  that  is  required  is  to 
retain  and  extract  this  information. 

Second,  integral  equation  synthesis  should  be  extended 
to  geometries  other  than  spherical.  Of  course,  as  it  is 
formulated,  the  method  is  independent  of  geometry.  However, 
the  transformation  of  the  theoretical  formulation  into  a 
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numerical  technique  has  been  accomplished  only  for  spherical 
geometry.  Potential  applications  of  this  efficient  technique 
will  be  greatly  enhanced  when  rectangular  and  cylindrical 
systems  can  also  be  treated. 
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Appendix  A 


Numerical  Methods  Used  in  the  IES  Calculations 


The  integral  form  of  the  neutron  transport  equation 
for  one  energy  group  can  be  written 


-/  E(r  +  pOdp- 

0  _ap 

*(r)  =  dV  ^ -  E(r ')C(r")^(r ')  e  v  (A.l) 

4irp  2 


where 


<K?) 

(?) 

(?0 

p 

Z(r') 

C(rO 


v 

a 


scalar  neutron  flux  in  neut/cm2-sec 

space  point  to  be  calculated 

other  space  points  over  which  we  integrate 


total  macroscopic  cross  section  at  r' 

average  number  of  neutrons  everging  from 
a  collision  at  ?' 

neutron  velocity 

neutron  growth  rate 


The  basic  problem  is  to  solve  Eqn  (A.l)  numerically 
and  efficiently  for  a  spherical  system.  There  are  many 
ways  in  which  this  can  be  done;  however,  the  following 
method  is  quite  efficient. 
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The  first  decision  involves  the  geometry  of  the 
integration.  So  far,  we  have  just  indicated  an  integration 
over  dV'  .  The  usual  procedure  is  to  define  the  vectors 
r  and  r'  from  the  center  of  spherical  geometry.  By 
appropriate  transformations  of  the  variables,  it  is  possible 
to  arrive  at  an  integral  equation  whose  kernel  is  symmetric 
in  r  and  r"  .  There  are  several  advantages  to  this 
approach,  particularly  in  that  the  solutions  to  the  equation 
then  form  an  orthogonal  set.  Also,  many  matrix  manipulation 
schemes  are  simpler  when  the  matrix  is  symmetric.  The 
difficulty,  however,  is  that  a  discontinuity  arises  when 
r'  *  r  .  This  can  be  avoided  by  sacrificing  symmetry  and 
using  p  as  the  variable  of  integration.  We  then  have 


dV  =  2 it  sinQdO  p2-  dp 


(A. 2) 


where  0  is  the  angle  between  r  and  p  ,  and  the 
discontinuity  is  eliminated. 

The  integral  over  dV'  now  involves  only  two  dimensions 
since  the  problem  is  symmetric  in  azimuthal  angle.  The 
integration  over  0  can  be  quickly  transformed  into  one 
over  p  =  cos©  since 


/  sin0d0  f(0) 


/  dp  f(p) 

-1 


(A. 3) 


t 


Integrals  of  this  type  are  most  accurately  approximated 
by  Gaussian  quadratures: 

1  K 

/  dp  f(p )  -  l  w,  f(p.)  (A. 4) 

-1  k=l  K  K 

For  each  value  of  K  ,  tables  (Ref  1)  give  the  corresponding 
values  of  the  weights  and  the  direction  cosines  p^ 

Using  Gaussian  quadratures,  the  entire  volume  integral  now 
reduces  to  a  weighted  sum  of  K  one-dimensional  integrals 
over  p  .  Integral  equation  synthesis  solutions  using 
K  =  10  have  been  shown  to  give  acceptable  convergence 
of  the  0-  integration. 

n  -jVcpOdp' 

K  p  o 

Ur)  =  l  W/^dpe  E(r')C(r')*(rO  (A.5) 

k=l  o 


where 


The  final  integration  over  p  requires  some  care.  We  must 
first  define  the  manner  in  which  known  properties  (total 
cross  section  and  average  number  of  neutrons  emerging  from 
a  collision)  and  unknown  quantities  (fluxes)  vary  as  functions 
of  radius.  We  must  then  establish  a  correspondence  between 
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this  radial  dependence  and  the  dependence  on  the  coordinate 
p  .  Because  each  integration  ray  is  oblique  to  the  radius 
vector,  this  latter  correspondence  can  become  geometrically 
and  numerically  quite  complicated. 

Before  describing  the  radial  dependence  of  the  variables, 
the  following  indices  are  introduced: 


k  =  the  label  of  the  angular  ray  along  which  we 

are  integrating  over  p  .  ; 

i  =  the  label  of  the  spherical  shells.  These  are 
labelled  such  that  i  increases  as  we  go  out 
to  larger  radii. 

j  =  the  label  of  the  crossing  points;  i.e.,  the 
jth  shell  intersected  as  we  proceed  outward 
along  the  kth  ray.  The  point  from  which  we 
start  is  j  =  0 

t  =  an  "offset"  index  that  is  a  function  of  k 
and  j  .  After  crossing  j  shells  as  we 
go  outward  on  the  kth  ray,  we  intersect 
the  (i  +  fc)  shell. 

i '  =  another  "offset"  index  which  is  either  the  same 

as  l  (if  the  kth  ray  is  going  inward  or  is 
just  "glancing"  across  a  shell)  or  equal  to 
(t  +  1)  if  the  kth  ray  is  going  outward.  This 
index  specifies  the  cross  sections  to  be  used 
in  the  integration. 


| 

i 


i 


An  example  of  these  indices  and  relationships  are  indicated 
in  Fig.  A-l. 


The  radial  dependence  of  the  various  properties  and 
quantities  is  approximated  as  follows: 

1.  Given  parameters  (cross  section  and  neutrons  per 
collision)  are  assumed  to  be  uniform  between  the  shell 
boundaries.  They  are  labelled  with  the  index  of  the  outer 
shell  boundary. 

2.  Fluxes  are  defined  oil  the  shell  boundaries  (that 
is,  at  specific  radii)  and  vary  linearly  with  r  between 
shell  boundaries.  We  see  that  the  actual  problem  is  repre¬ 
sented  by  an  approximate  configuration  in  which  the  neutronic 
parameters  are  uniform  in  concentric  annuli,  whereas  the 
fluxes  vary  linearly  with  r  within  a  given  shell.  We 
would  expect  this  approximation  to  be  valid  insofar  as  the 
mean  free  path  is  long  compared  to  the  thickness  of  the 
shell.  Thus,  the  larger  the  mean  free  path  relative  to  the 
scale  of  the  geometry,  the  more  accurate  is  this  representation. 

The  integration  along  each  ray  requires  that  the  flux 
be  expressed  as  a  function  of  p  rather  than  r  .  It  is 
convenient  to  represent  1K0)  by  means  of  interpolation 
coefficients ,  cjj^  (Ref  17).  The  upper  indices  indicate 
the  position  in  the  mesh:  we  started  at  radius  i  and 
have  progressed  outward  on  ray  k  ;  this  interpolation  will 
apply  between  boundary  crossings  j  and  j  +  1  .  The 

lower  indices  relate  to  a  polynomial  fit  involving  the 
values  of  flux  at  radii  (i  +  l  +  m)  ,  the  powers  of  (■—-) 
being  given  by  n  .  These  terms  are  illustrated  in  Fig.  A- 2. 
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For  p  between  pik^  and  Pikj+1  ,  the  flux  is 
given  by 

=  -1,0  ^i+5,-1  +  -1,1  ^i+ji-l  ^  +  -1,2  ^i+fc-l  ^ 

+  C0,0  ^i+i,  +  C0,l  0i+«.  ^  +  C0,2  *i+£  ^ 


+  Cl,0  ^i+£+l  +  Cl,l  ^i+t+l  Q  +  Cl,2  ^i+2,+1  ^  (A*6) 


or 


Mp) 


l 


I 

m=l 


2 

V  C  <i. 
n=o  m,n  1+i+m 


V 


n 


(A.  7) 


where  the  upper  indices  have  been  suppressed  for  the  sake 
of  clarity. 

The  calculation  of  these  interpolation  coefficients 
is  a  somewhat  tedious  algebra  problem,  but  once  computed 
they  can  be  stored  and  re-used  for  a  given  radial  mesh. 

A  considerable  simplification  is  possible  if  the  radial 
mesh  is  made  self-similar  wherein  each  radius  r^  is  a 
constant  factor  (e  >  1)  times  the  preceding  radius  r^  ^ 
In  this  case,  the  coefficients  C  are  no  longer  dependent 
on  i  ,  the  label  of  the  initial  radius. 

The  interpolation  coefficients  are  calculated  from 
Eqn  (A. 7)  with  the  following  boundary  conditions: 
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@  p  =  Pj 


(A. 8) 


;  *(p)  =  *i+* 


§ 


djp_  s_r 
8r  op 


(A. 9) 


@  p  1=1  Pj+1  ; 


<p 


i  +  Z+l 


(A. 10) 


Note  that,  depending  on  the  angle  of  the  integration  ray  k  , 
Pj+^  will  be  on  a  shell  whose  radius  is  larger  than  (outward 
case),  the  same  as  (glancing  case),  or  smaller  than  (inward 
case)  the  shell  denoted  by  p^  .  For  this  reason,  the 
interpolation  coefficients  are  also  divided  into  outward, 
glancing,  and  inward  cases.  A  listing  of  the  interpolation 
coefficients  for  a  self-similar  radial  mesh  is  provided 
in  Table  A-l. 

At  this  point,  making  use  of  Gaussian  quadrature  and 
the  interpolation  coefficients,  we  can  rewrite  Eqn  (A.l)  as 


•'max 


’H1’)  I  I  I  1  Cjjjn  '^+£+m  ^i+#, 

k=l  K  j=0  m=-l  n=0  m  1  *  m  14  1  £ 


r 


Pj+1 

/  dp  (£-) 

P-i  1 


V 


J 


(A.  11) 
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The  terms  above  in  the  brackets  are  simplified  in 

-ap/v  „„  M  _  ap 


the  following  way.  Approximate  e 
Define  a  new  variable  of  integration  as 


as  Cl  -  • 


p  -  p. 


y  - 


pj+i  ‘  pj 


i  Ic 

Define  a  transmission  factor  T.’  such  that 

J 


rpi  ,  k 

j+1 


T1,k  e 


-r 


j+1  ^pj+l  ‘  pj 


and 


■pi  >k  „ 
0 


Then  the  bracketed  term  in  Eqn  (A. 11)  can  be  written  as 


=  T 


J,k  /  dy  yn  [1  -  £  Upy  -  p.)]  e  1+)l  (A.l 


2) 


where 


Ap  E  pj+l  '  pj 


Table  A-l 

Interpolation  Coefficients  for  A  Self -Similar  Radial  Mesh 


M  Case 


OUTWARD 


GLANCING 


INWARD 


n  K  -  XjXj+l 


X  -X  .  X  . 

-  1  JL+i+_i. 


^  ‘  gmin^ 


n  2  -  x  •  x .  . 


AX  AX' 


OUTWARD  I  1+g " 


X  -X  .  ,  X . 

11+1  1 


Ax  Ax' 


0  GLANCING  K 


INWARD 


.  'VxJ*i>. 

2xj 

a 

Ax 

AX2 

(x . +X .  , ) 

I  Jt.l  X 

2Xj 

5 

AX 

AX2 

■8  - 2  8 

AX  AX2 


AX  AX' 


AX  AX' 


X  .  X  «  x  . 

OUTWARD  -g'.  1  i-+-+  J— 
AX  AX2 


fx.+x.A,)  2x.  ..  , 

_j _ mi .  —±  -g'A.  +  j. 

AX  AX2  AX  AX' 


GLANCING 


Table  A-l  (Cont'd) 


=  cose  where  0  is  the  angle  between 
and  the  kth  ray 


x . 


J 


'K 


Then  the  Table  entries  are 


rj 

m,n 


59 


A  second  change  of  variables  such  that 


z  -  ap  y 


allows  us  to  write  the  bracketed  term  as 


(1  +  -  o.)  (  -  )n+1  /  Zn  e_Z  dZ  -  (  —  Ap)  C  —  )n+2 

v  J  b  o  v  b 


/  Zn+1  e'Z  dZ 


(A. 13) 


where 


b  =  Z.  +  Jl  Ap 


But  these  integrals  are  just  the  partial  gamma  function 
such  that 


/  e"Z  dZ  =  e"a  e"b  =  G  (a,b) 
a  o 


/  Zn  e’Z  dZ  =  -bn  e"b  +  n  /V'1  e'Z  dZ 


CA. 14) 
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or 


Gn(a,b) 


'n  -Z  . 

Z  e  dz 


=  -b 


n 


n  Gn-lfa,b^ 


Making  these  substitutions  and  rearranging,  we  have  the 
final  numerical  form  of  Eqn  (A.l): 


*i 


K  Jmax  1 

l  l  l 

k=l  j=0  m=-l 


•pi  ,k 

'  j 


rj,k 

Lm,n 


Gn(o,b) 


v. 


Before  this  equation  can  be  solved  numerically,  the 
boundary  conditions  need  to  be  incorporated.  Incorporating 
the  external  boundary  condition  is  a  trivial  matter;  the 
summing  over  j  (the  number  of  shell  crossings)  continues 
until  the  outermost  boundary  is  reached.  A  simple  test  of 
(i  +  i  +  m)  denotes  when  the  outer  boundary  is  reached. 

The  interior  boundary  condition  is  a  little  more  complicated. 
Because  of  spherical  symmetry,  the  internal  boundary  condition 
is 


(A. 16) 
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Obviously,  a  linear  variation  of  flux  with  radius  inside 

the  smallest  shell  (r  <  r  )  cannot  satisfy  the  interior 

boundary  condition.  For  this  reason,  a  quadratic  variation 

of  flux  with  radius  is  assumed  inside  r  .  For  r  less 

o 

than  r 

o 

Hr)  =  A  +  Br  +  Cr2  (A. 17) 

Denoting  the  flux  at  the  center  Hr  =  0)  as  'J'q  and 
applying  Eqn  (A. 16),  we  have 

ip  -  ip 

-Hr  <  r  )  =  H  -  -Z - £  r2  (A. 18) 

°  L  r  2 

O 

where  ip  =  Hr„)  .  The  center  flux  i>r  is  computed 

O  O  L * 

implicitly  as  a  function  of  all  the  ip^  .  That  is, 

I 

max 

1>r  ■  I  a.  (A. 19) 

L  i=0  1  1 

where  the  a^  are  obtained  from  a  solution  of  Eqn  (A.l) 
at  r  =  0  .  Even  though  Eqn  (A.l)  is  greatly  simplified 
when  r  =  0  ,  the  resulting  expansion  coefficients  a^^ 

are  algebraically  complicated.  They  are  presented  here 
without  derivation. 


C  G-(e  r  r  )  ,  C, 

ao  =  TT - ~ - 2~2"  4  E  TF~-r  T  [rlGo(£lro>Elri)Tl 

r  i  1  o '  °  1  1  1 


C.  T. 


ai  = 


TrT-TT — 5 - IT”  (ziri_i  >  Siri^  "  ri-i^o  ^iri-l  ’ 


T.  C . 

*  DTr^-r,)  ^ri+lGo^£i+lri,£i+lri+l^ 


Gl(Ii+lri'Ei+lri+l 


)3 


(A. 20) 


where 

D  ’  [r-  -  W-i’W  *  r^W-i-Vi^o 

ro 

Tj  ■  eJtPCri-l(Ji-  =k-l)]Ti.l 

To  ’  1 

6n  *  partial  gamma  function  as  defined  in  Eqn  (A. 14) 


Having  incorporated  the  boundary  conditions,  the 
transformation  of  Eqn  (A.l)  to  a  discretized  form  is  now 
complete.  Clearly,  Eqn  (A. 15)  is  an  eigenvalue  problem  of 
the  form 
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(A. 21) 


It  can  be  solved  numerically  by  techniques  such  as  the 
power  method  (Ref  7:291).  The  procedure  is  to  solve  Eqn 
(A. 15)  and  then  adjust  a  and  resolve  the  equation  until 
the  eigenvalue  is  made  equal  to  unity.  The  a  that  makes 
the  system  just  critical  is  the  one  group  growth  rate  that 
we  seek. 


64 


Appendix  B 


Data  Listing  for  Jezebel  Calculations 

Jezebel  calculations  were  performed  with  the  IES 
method  using  a  17  shell  self-similar  mesh  with  rQ  =  0.5354 
and  e  =  1.157129  .  The  Gaussian  quadrature  abscissas  and 

weights  are  listed  in  Table  B-I. 


Table  B-I 

Gaussian  Quadrature  Abscissas  and  Weights3 

Abscissas 

0.97391 

0.86506 

0.67941 

0.43340 

0.14887 

Weights 

0.06667 

0.14945 

0.21908 

0.26927 

0.29552 

Abscissas 

-.14887 

-.43340 

-.67941 

-.86506 

-.97391 

Weights 

.29552 

0.26927 

0.21908 

0.14945 

0.06667 

3  Data  from  Ref  1:916. 

The  sixteen  energy  groups  and  their  associated  average 
group  velocity  and  cross  section  used  in  Jezebel  calculations 
are  provided  in  Table  B-II. 
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Table  B-II 


Group  Structure,  Velocities  and  Cross  Sections 
For  Jezebel  Calculation3 


Group 


Group  Energy 
Lower  Bound  (MeV) 

Group  Velocity 
(cm/shake) 

Transport  Cross 
Section  (Barns) 

13. 62b 

51.60 

3.8553 

12.57 

50.24 

3.8490 

11.57 

48.19 

3.6706 

9.68 

45.24 

3.5983 

7.97 

41.00 

3.8157 

6.42 

37.01 

4,0943 

4.73 

32.48 

4.2653 

2.55 

25.68 

4.7049 

1.18 

18.26 

5.0878 

0.638 

12.91 

5.8035 

0.296 

9.263 

6.8698 

0.  222 

6.988 

8.4536 

0.107 

5.373 

9.7588 

0.0273 

3.  229 

11.4430 

0.00875 

1.739 

13.9460 

0.00100 

0.7801 

16.5620 

From  Ref  22. 

Upper  Energy  bound  on  group  1  is  14.60  MeV. 
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